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Abstract 

Irregular bone remodeling is associated with a number of bone diseases such as osteoporosis and 
multiple myeloma. Computational and mathematical modeling can aid in therapy and treatment as 
well as understanding fundamental biology. Different approaches to modeling give insight into different 
aspects of a phenomena so it is useful to have an arsenal of various computational and mathematical 
models. Here we develop a mathematical representation of bone remodeling that can effectively describe 
many aspects of the complicated geometries and spatial behavior observed. 

There is a sharp interface between bone and marrow regions. Also the surface of bone moves in and 
out, i.e. in the normal direction, due to remodeling. Based on these observations we employ the use of a 
level-set function to represent the spatial behavior of remodeling. We elaborate on a temporal model for 
osteoclast and osteoblast population dynamics to determine the change in bone mass which influences 
how the interface between bone and marrow changes. 

We exhibit simulations based on our computational model that show the motion of the interface 
between bone and marrow as a consequence of bone remodeling. The simulations show that it is possible 
to capture spatial behavior of bone remodeling in complicated geometries as they occur in vitro and in 
vivo. 

By employing the level set approach it is possible to develop computational and mathematical repre- 
sentations of the spatial behavior of bone remodeling. By including in this formalism further details, such 
as more complex cytokine interactions and accurate parameter values, it is possible to obtain simulations 
of phenomena related to bone remodeling with spatial behavior much as in vitro and in vivo. This makes 
it possible to perform in silica experiments more closely resembling experimental observations. 
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Introduction 



Bones of the mature skeleton undergo constant remodeling through resorption and replacement of old matrix. 
The repair of micro-cracks caused by normal mechanical stresses entails stimulation of local remodeling by 
tissue-intrinsic factors [3,19,25]. Bone remodeling also plays a role in systemic Ca ++ homoeostasis and is 
globally regulated by humoral factors (e.g. parathyroid hormone). Because the major mechanisms regulating 
bone metabolism are well-established and the links between abnormal metabolism and bone diseases are 
well-documented, the remodeling process is a potentially fertile ground for mathematical and computational 
analyses. Computer models that simulate the complexly regulated biology of bone remodeling can be used to 
rapidly and inexpensively screen drugs to estimate their disease-modifying activity. By minimizing the need 
for much more lengthy and expensive in vivo testing, such models can greatly accelerate the development of 
novel treatments for bone diseases. 

Bone remodeling occurs simultaneously at various locations. However, remodeling from one site to another 
is not necessarily synchronous. Remodeling is localized to a particular remodeling site with teams of cells 
acting together as an individual structure commonly referred to as a basic multicellular unit (BMU) [3,19,24, 
25]. Remodeling is driven locally by signals from osteocytes, which secrete chemotactic factors in response 
to mechanically induced microfractures and in response to bone strain [3, 19,25]. With the release of factors 
by osteocytes stromal stem cells are activated and begin to produce factors, such as macrophage colony- 
stimulating factor (M-CSF), that promote osteoclastogenesis [3,19]. Additional stromal cells are created and 
differentiate into pre-osteoblasts. These cells begin producing factors such as receptor activator of nuclear 
factor-B ligand (RANKL), which binds to receptors present in pre-osteoclasts and inhibits apoptosis [3,19]. 
This is a case of osteoblast-derived paracrine signaling and is featured in many of the mathematical treatments 
of remodeling discussed below. In the presence of the factors M-CSF and RANKL pre-osteoclasts fuse to 
become mature osteoclasts. Osteoclasts bind to bone surfaces and resorb the matrix by secreting acids 
and matrix proteases. Mature osteoclasts in turn produce growth factors and other signaling molecules, an 
example of osteoclast-derived autocrine signaling. Resorption also releases growth factors such as IGF and 
TGF-/3 that act on cells of the osteoblast lineage, an example of osteoclast-derived paracrine signaling [3,19]. 

When resorption is nearly complete pre-osteoblasts mature into osteoblasts and stop secreting RANKL. 
Instead they produce OPG, a decoy receptor that blocks RANKL binding to osteoclasts. As a result, 
osteoclasts cease resorptive activity and undergo apoptosis. During formation some osteoblasts become 
trapped in the matrix as osteocytes, while some undergo apoptosis [3, 19]. This completes the description of 
a typical scenario for remodeling by a BMU in response to mechanical stress. This description is helpful in 
the development of the mathematical and computational models discussed below. 

A first step towards formulating a theoretical representation of bone remodeling is to incorporate the 
mechanisms that regulate the process. Here this is the autocrine and paracrine biochemical signaling amongst 
cells. Komarova et al [11,12] develop a temporal model consisting of ordinary differential equations (ODEs) 
describing the local dynamics of the remodeling process. This model uses a power law approximation to 
capture the interaction of cell populations through autocrine and paracrine signaling. An alternative approach 
to that in [11, 12] also capturing the cytokine signaling regulating bone remodeling is given by Pivonka et al 
in [23]. The model presented in [23] incorporates signaling through Hill functions rather than power laws. 
Both approaches model the changes in osteoclast and osteoblast populations in response to RANKL/OPG 
signaling. Ayati et al [2] and Ryser et al [26, 27] build on the model from [12] for the local dynamics to 
incorporate some of the spatial effects of bone remodeling. In Ayati et al [2], in the context of multiple 
myeloma, diffusion of cells is added to incorporate some spatial effects of bone remodeling, while in [26] 
Ryser et al consider chemotaxis of osteoclasts in the presence of remodeling promoting factors. 

In this article we employ a geometric method developed by Oshcr and Sethian [17] to present another 
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approach to capturing spatial effects of bone remodeling. We represent, via the so called level-set equation, 
the evolution of the interface between bone and marrow regions in response to remodeling. The interface is 
described geometrically as a closed curve in the plane. The interior of the curve models region of bone while 
the exterior to the curve represents the marrow region. A level-set function, which is determined by solving 
the level-set equation, provides a convenient way to represent a geometric object such as a curve. Through 
this geometric approach it is possible to consider effects of remodeling for complex geometries and spatial 
scales not easily represented with previous spatial models. This is particularly relevant for representing the 
remodeling of cancellous or trabecular bone tissue. By interpreting the change in total bone mass, at each 
point in space, during remodeling as the speed of the interface motion (normal, i.e. perpendicular, to the 
interface) 1 we obtain a novel representation of the spatial behavior of bone remodeling. We note that while in 
many other applications of this geometric approach, the motion of the interface is determined by some aspect 
of the local geometry such as curvature, this is not the case for the application we give here. The motion 
of the interface is a result of the behavior of the cell populations at each point in space, in other words, the 
changes in geometry are directly determined by biological phenomena alone. We find that our modeling and 
simulation efforts are representative of observations in animal models. Incorporating pathological behavior 
associated with abnormal bone remodeling in future work is straightforward and of clinical interest. 



Models 

The complex nature of the interactions of osteoclasts and osteoblasts and the chemical signaling that regulates 
them makes a detailed theoretical and experimental study difficult. Mathematical modeling may be used 
to tie existing knowledge about a system with new ideas. Simulations based on mathematical models can 
provide a useful description of observed or hypothesized results otherwise difficult with ordinary language. 
In this section we present and discuss some mathematical models of bone remodeling with a focus on the 
ones that have influenced the development of our approach. The models presented here describe the change 
in osteoclast and osteoblast populations throughout the remodeling process. Throughout we denote by C 
the osteoclast population and B the osteoblast population. 

The model which forms the foundation for all that follows is the temporal model of Komarova et al [12] 
that describes the local dynamics of bone remodeling at a single remodeling site. The dynamics of the cell 
populations is determined by the interaction of cells through autocrine and paracrine signaling. In [12] the 
authors employ a power law approximation to describe the overall effects of the biochemical factors involved 
in signaling. This type of modeling, while simple to describe, is effective in capturing the nonlinear behavior 
of biochemical signaling. The exponents in the power law approximation describe the effectiveness of the 
autocrine and paracrine signaling on each cell type. From [12] the dynamics of the cell populations are given 
by 

dC 

— = a 1 C»"fi» I -0iC, (la) 

— = a 2 C Sl2 B 922 - (3 2 B. (lb) 
at 

Here a\ and a 2 are constants that reflect the rate of cell production while while (3i and (3 2 are constants 
reflecting the rate of cell removal of osteoclasts and osteoblasts respectively. As mentioned above the expo- 
nents gij describe the effectiveness of the different types of signaling on each cell type with gn describing the 

iMore precisely, at each point of the interface one can put a tangent line. By normal we mean perpendicular to this tangent 
line. 
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effectiveness of osteoclast derived autocrine signaling, 1721 the effectiveness of osteoblast derived paracrine 
signaling, 1712 the effectiveness of osteoclast derived paracrine signaling and 1722 the effectiveness of osteoblast 
derived autocrine signaling. The authors use this model to simulate the behavior of remodeling and find that 
the value of the parameter gn for the effectiveness of osteoclast derived autocrine signaling most strongly 
effects the mode of remodeling behavior. 

The authors also calculate the change in total bone mass due to remodeling via the equation 

dz 

— = -k x max{0, C - C] + k 2 max{0, B - B}. (lc) 
at 

Here z is the bone mass k\ , hi are constants describing the activity of bone resorption and formation 
respectively. The constants C and B denote the stable steady state values of osteoclasts and osteoblasts 
respectively and here are given by 

1-922 „ „ 921 

\aij V a 2/ 

(. 912. , . 1-911 

M ' (A) ' , 

7 = .912.921 - (1 - .911 )(1 - .922) 

In the level-set representation that follows we take the change in bone mass as giving the speed of motion 
of the interface in the normal direction. We note that in the work of Pivonka et al [23] there is an equation 
for the change in bone volume similar to (lc). Thus the equations from [23] can be used in the following in 
place of (la) and (lb). However, for computational simplicity we use the power law approximation since the 
resulting interface dynamics should be the same, for the same z field, regardless which model is used for the 
cell populations is used to obtain such a z field. For example, in our case we can obtain the same z field 
from a system with and without explicit chemical concentrations. For convenience we reproduce (see figure 
1) the change in bone mass computed in [12] as a percentage of initial bone mass using the system (la), (lb), 
and (lc). Figure 1 can then be used for qualitatively comparison with the results in this paper. 

In [2] the authors give an implicit spatial representation of bone remodeling by interpreting the bone mass 
equation as reflecting bone thickness. We use this interpretation with the level set method discussed below 
to give a first study of the dynamic behavior of the bone marrow interface due to remodeling. This assumes 
an unrealistic degree of homogeneity but is an instructive first step in developing the level set approach. The 
authors in [2] go on to incorporate explicitly spatial effects through the addition of linear diffusion terms 
to the system (1). This is all done in the context of multiple myeloma bone disease since the goal of their 
work was to show how perturbation of the mathematical models for bone remodeling can exhibit pathological 
behavior as is seen in association with metabolic bone diseases. The modeling approach developed below 
can also be used to study how pathological behavior can arise from perturbations representing irregular 
remodeling. 

In [26] the authors develop a spatial representation of bone remodeling by modifying the system (1) to 
include motility of osteoclasts and explicit representation of the RANKL and OPG fields. The model is 
particularly successful in capturing the motion and steering of the BMU and simulations in one and two 
dimensions show the distinct features of the cutting cone of osteoclast cells. Their results are obtained by 
including the chemotaxis of osteoclasts towards increasing concentrations of RANKL and explicitly the dy- 
namics of the RANKL/OPG pathway. The equations for RANKL and OPG concentrations include porous 
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Figure 1: Change in bone mass according to the model from [12]. This figure was obtained by solving 
equations (la), (lb), and (lc). Here the horizontal axis is time in days and the vertical is percentage of initial 
(100%) bone mass. 

flow of the chemicals, activation of osteclastogencsis by RANKL and binding of RANKL/OPG. The two di- 
mensional simulations based on the model elegantly show steering, as a chemotactic response, of a remodeling 
BMU along a microfracture . 

Also of note are the mathematical models of Martin and Buckland- Wright [13,14]. The model in [13] 
simulates erosion depth due to resorption by modeling the changes in cellular activity due to cytokine 
signaling. This is done using the equations for enzyme kinetics of Michaelis and Menten. The model in [13] 
differs from the mathematical models in [2, 11, 12,23,26] in that it does not describe changes cell populations 
while the model in [14] does include changes in the osteoblast population. The models found in [13, 14] 
provide another useful method for studying how changes in signaling are manifested as physical changes 
in the bone. The recent review paper of Pivonka and Komarova [22] gives a nice overview mathematical 
modeling in bone biology. 

In the models discussed so far bone remodeling is regulated by biochemical signaling alone. Mechanical 
signaling also plays an important role in the regulation of cells involved in bone biology. There is a substantial 
body of work devoted to combining the biochemical and mechanical stimuli in mathematical models to study 
the effects that the two different types of stimuli have together [4-8,20]. 

For example in [8] the authors couple mechanical stimuli through finite element analysis with partial 
differential equation models describing spatio-temporal changes in cell populations. These changes are due 
to cell migration (chemotaxis, haptotaxis), proliferation and differentiation. They also include equations 
describing angiogenesis. Biochemical signaling is incorporated by including interactions of cell populations 
with growth factors. This all leads to a mathematical framework that is successful in simulating the combined 
effects of biochemical and mechanical regulatory mechanisms. An interesting result is the prediction of over- 
load induced nonunion formation. The modeling approach just discussed is well suited to the study of 
wound and fracture healing and issues such as angiogenesis that arise under such situations. This is the 
case for example with larger fractures in or on cortical bone, rather than microfractures. Here we leave the 
role of angiogenesis to future considerations and focus on the spatial manifestations of other features of the 
remodeling process. 
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The geometric approach mentioned, was initially developed by Osher and Sethian in [17] to study geometry 
driven motion of an interface and is particularly well suited for problems in which there is a definite "inside" 
and "outside" such as in the modeling of soap bubbles. The method has since been further developed and 
extended, see for example [16,18,29,31,32], and has found many applications in various areas of applied 
mathematics and engineering [16,29,32]. Since the surface of bone forms a distinct, sharp interface the level- 
set function seems well suited for representing physical changes to the surface of bone due to remodeling. 

We present here some basic background on the level-set function approach specific to representing the 
spatial effects of bone remodeling in trabecular bone tissue. For a more in depth treatment and other 
applications see the basic texts by Osher and Fcdkiw [16] and Sethian [29]. As stated previously, we would 
like to describe the interface between bone and marrow regions as a closed 2 curve in the plane. The basic 
idea is to represent an interface, geometrically a curve, surface, etc. implicitly as some level set, typically 
the zero level set, of a higher dimensional object described by a so called level-set function. For example 
consider the function cf>(x, y) := 1 — x 2 — y 2 . The zero level set of <p(x, y) is the set of all points (x, y) in the 
plane satisfying 

4>(x, y) := 1 - x 2 - y 2 = 0, (3) 

which describes a circle of radius 1 centered at the point (0,0). Alternatively one could represent an in- 
terface locally by explicit functions, but this may require many different such function, or parametrically. 
However, the approach to representing geometry via a level-set function provides some advantages over these 
alternatives that make it well suited for applications. 

One advantage to representing an interface, in this case a curve, implicitly as above is that most of the 
geometric information one would like to have about the interface can be obtained from the level set function. 
For more details on this see the first chapter of the book by Osher and Fedkiw [16]. What really makes the 
level-set function appealing for applications such as the changes to bone surface due to remodeling, is that, 
given an initial interface represented by a specified level set function, there is a simple way to evolve the 
interface, by evolving the level set function in time. More precisely, to represent the evolution of an interface 
we use a family of level-set functions <f>(x, t) parameterized by time t so that at each time t the zero level set 
{x. : (f>(x, t) = 0} represents the interface at time t. Here x £ M. n for some positive integer n, the case we are 
concerned with has n = 2. Then following Osher and Sethian in [17] we can derive a differential equation 
for the motion of the interface in the normal direction due to normal velocity with speed given by a function 
a{x,t). The equation is 

4> t (x,t) = a{x,t)\V<j){x,t)\ (4) 

and is known as the level-set equation. Here (f>t denotes the derivative with respect to time, \/(j)(x,t) is the 
gradient of </> in the spatial variables and |y| is the length of a vector y € M. n . In a sense the function a(x, t) 
gives the rule telling what drives the motion of the interface and solving the level set equation describes 
what the interface looks like at any given time. Notice that here the level set equation does not contain a 
curvature term so that local geometry is not used to move the interface. 

The ODEs (la) and (lb) effectively describe the local population dynamics of the primary cells involved 
in remodeling. Alternatively, these ODEs can represent a larger, but well-mixed spatial domain. We have 
discussed previous approaches to incorporating spatial heterogeneity into mathematical descriptions of bone 
remodeling. Even with the success of these models it is valuable to move towards new and different spatial 
representations of bone remodeling. In that direction we employ the level-set function approach to develop 
such a representation. The interface between bone and marrow is distinct and sharp and we use this to 

2 We note that the curve need not be closed to apply the level set method. We take a closed curve and in particular a circle 
to make clear what is "inside" and what is "outside" . 
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motivate representing the system as one with a moving interface. We use the level-set equation discussed 
previously to describe the motion of this interface. 

Since the power law approximation of [12] successfully captures the local dynamics of osteclast /osteoblast 
cell populations we elaborate (la) (lb) by adding explicit space dependence and coupling with the level-set 
equation to represent spatial manifestations of the remodeling process. In other words, we (locally) couple 
the biological aspects, i.e. population dynamics of remodeling cells, of bone bone remodeling with the change 
in geometry of the interface between bone and marrow regions. This is done by assuming that the speed in 
the direction normal to the interface, a(x, t) in (4), is given by the change in bone mass. More precisely, the 
speed in the normal direction is given by a constant multiple of the change in bone mass, 

a(x, t) = v{-ki max{0, C(x, t) - C} + k 2 max{0, B(x, t) - B}), (5) 

where v is a constant. We note that here there is a (potentially) different population of each cell type at each 
location x in space. This makes the speed a(x, t) different at each point in space. In [26] chcmotaxis plays the 
role of a steering mechanism, determining the spatial behavior of the osteoclast and osteoblast populations 
along an existing microfracture on the order of a few microns. In our model the spatial behavior, including 
steering of cells, is determined by the geometry via the level set method. This is appropriate for spatial 
scales on the order of a thousand microns or larger and is well suited for geometrical configurations common 
in trabecular bone. If necessary, slight modification of our approach allows for the inclusion of chemotaxis 
or other steering mechanisms. 

Results and Discussion 

In what follows, we denote by T t = {x : 4>(x, t) = 0} the interface at time t, where is the level-set function. 
We solve the following system of equations to simulate the motion of the interface between bone and marrow 



due to bone remodeling, 
dC 

— (x,t) = a 1 C(x,t) 311 B(x,t) 921 -p x C(x,t) (6a) 
dB 

— (x,t)=a 2 C(x,t) 9l2 B(x,t) 9 ™ -p 2 B(x,t), (6b) 

(j)t = a(x,t)\V(t>\, (6c) 
C(x,t) = C(x,t), xeV t (6d) 

B(x,t) = B(x,t), x€T t (6e) 



where a(x, t) is given by (5), and C(x, t) and B(x, t) are the number of osteoclasts and osteoblasts respectively, 
recruited from the marrow to a point x on interface at time t. Thus we have a system of equations where 
C{x,t), B(x,t) depend on the interface, i.e. level set function (j)(x,t), since at each time, the population at 
the interface T t is updated according to (6d) and (6c). The rapid motion of the cells from deeper within the 
marrow to a remodeling site is captured by the local recruitment terms (6d) and (6c). When remodeling is 
initiated the cells are recruited quickly and directly to the remodeling site. As a result, diffusion of these 
cells is not relevant for the purposes we are considering. The domain for this problem is a square. To mimic 
a large domain, for computational simplicity we assume periodic boundary conditions for this domain. The 
initial interface is described by the zero level set of an initial level-set function which we denote by 4>o(x,y). 
For this discussion we take the initial interface to be a circle but note that any other shape could be taken by 
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prescribing an appropriate level-set function. It is even possible to use, as a level set function, images taken 
from experiments, see for example [30]. For the numerical simulations presented below we use parameter 
values as in [2, 12,27], these are 



C = 1.06 


B 


= 212.13 


ai = 3.0 


Pi 


= 0.2 


a 2 = 4.0 


$2 


= .02 


= 0.5 


921 


= -0.5 


.912 = 1.0 


922 


= 0.0 


fci = 0.24 


fc 2 


= 0.0017 



All of our simulations have been carried out using MATLAB. The implementation of our model required 
solving a coupled system of differential equations at each point in the domain. To carry this out we first 
solved the equations for the cell populations using implicit time stepping and a Krylov iterative method [9,10] 
to solve the nonlinear system. This provided the normal speed for the motion of the interface. Evolving the 
interface required solving the level-set equation. For this we used the level set toolbox for MATLAB [15] 
developed by Ian Mitchell. This toolbox contains many useful tools for implementing the level-set numerical 
method in MATLAB. 

We begin remodeling at three distinct sites by perturbing the osteoclast population away from steady 
state, figure 2(a). This initiates bone resorption. The effects of remodeling on the interface after 5 days of 
bone resorption is shown in figure 4(b). Note that the interface moves inward (toward the interior region 
on bone) due resorption since there is an increase in osteoclast population away from steady state. The 
osteoclast population quickly returns to steady state, figure 2(f). In accordance with equation lc when the 
osteoclast population returns to steady state the motion should cease to be inward. This can be seen in 
figures 4(e) and 4(f). In the meantime the osteoblast population increases from steady state, figures 3(b), 
3(c), and new bone is being formed. The effect of remodeling on the interface during the early stages of 
formation of new bone by osteoblasts is shown if figures 4(f) and 4(g) where the indentions become filled 
in. After reaching approximately maximum population density, figures 3(e), 3(f), the osteoblasts return to 
steady state by which time bone has been fully remodeled, figures 3(i), 4(i). 

Comparing the results in the three figures 2,3 and 4 at corresponding times shows how the changes in the 
cell populations are manifested as physical changes in the bone surfaces at three different remodeling sites. 
The interface moves in the normal direction as a result of remodeling in accord with (1). While here we have 
initiated remodeling at each of the three sites at the same instant this is not a requirement of this modeling 
approach. Remodeling at different sites can easily be taken to be out of phase. 

We have described in figure 5 quantitatively the change in the bone by plotting the area (in square 
microns) interior to bone, this can be compared with figure 2, which shows the bone mass as a function of 
time. We remark that bone mass resorption and formation occurs somewhat faster in figure 5 as we add 
more spatial dimensions due to the geometry, (i.e. proportional to radius squared versus proportional to 
local thickness). We note however, the qualitative agreement with figure 2 in [12]. 

As discussed previously, in the work [26], a mathematical model including explicit spatial terms is con- 
sidered. To illustrate the flexibility of the level-set approach we discuss how this can be done here as well. 
The first step is the inclusion of signaling pathways not through exponents in a power law approximation 
but explicitly. For example the RANKL/OPG pathway can be represented explicitly by modifying (la) to 
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Figure 2: The population dynamics of osteoclast cells at three separate remodeling sites. Re- 
modeling is initiated at three distinct sites by perturbing the osteoclast population away from steady state. 
The osteoclast population quickly returns to steady state. Bone is resorbed in the presence of an osteoclast 
population above steady state levels. 

include the effect of RANKL on osteoclast arriving at 



t = aiC 911 — fiiC + Ki . 
at A + u 



H (max{0, C - C})C. 



(7a) 



Then the RANKL/OPG field is represented by equations for the concentrations of each chemical species. 
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Figure 3: The population dynamics of osteoblast cells at three separate remodeling site. The 

osteoblast population starts at steady state, then increases in response to an increase in osteoclast population, 
and finally returns to steady in accord with (la) and (lb). Bone is formed in the presence of an osteoclast 
population above steady state levels. 



These are 



u t = D(x)Au - K 2 - (max{0, C - C})C 

A + u 

+ an max{0, B — B} — n^uv, 
v t = D(x)Av + ao max{0, B — B} — k^uv, 



(7b) 
(7c) 
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Figure 4: The evolution of the interface between bone and marrow regions due to remodeling 
at three distinct remodeling cites. The speed in the normal direction is given by the change in the bone 
mass (5). Figure 4(a) shows the initial interface before any remodeling has begun. The interior of the circle 
represents the bone region while the exterior is marrow. The interface between bone and marrow moves in 
and out, i.e. in the normal direction, due to changes in the populations of osteoclasts and osteoblasts. 



where u denotes the concentration of RANKL and v the concentration of OPG. The equations for the 
osteoblast population and change in bone mass remain as in (lb) and (5) respectively. Note that while 
in [26] steering is controlled by the chemotaxis of osteoclasts toward RANKL concentrations, here it is 
controlled by the geometry through the level-set function as above. Our computational results including 
explicit spatial terms for the chemical species showed an evolution of the interface similar to figure 4. 
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Figure 5: The amount of change due to remodeling. The units are days for the horizontal axis and 
mm 2 for the vertical axis. We arrived at this by computing the interior area at each time in the evolution 
of the interface. This can be compared with figure 1. The differences in geometry (i.e. velocity times radius 
squared here, and local bone mass in 1) results in a slightly faster, but qualitatively similar, cycle than in 
figure 1. 

Based on the local dynamics of bone remodeling as developed in [11, 12, 26, 27] we have developed a new 
approach to modeling the spatial effects of bone remodeling. Through the use of the level-set equation we 
link local dynamics with geometric effects and spatial behavior such as steering of a BMU in a way that 
is appropriate on many different spatial scales. The approach taken is useful but not limited to describing 
the effects of bone remodeling in complex geometries such as is seen in the remodeling of trabecular bone. 
The extension of this model to various other situations is straightforward and software for implementing the 
level-set equation approach numerically is widely available (see e.g. [15] and references therein). 

The development of mathematical and computational models of biological phenomena is valuable in 
design and analysis of experiments and aids in building intuition and conceptual understanding of biological 
processes. These models can also be useful to the biomedical community. Analysis and simulation based on 
mathematical and computational models can suggest how and where pathological behavior occurs. Examples 
of this in bone remodeling can be seen in [2, 11, 12] and the modeling approach in this paper can be used in 
the same way by including the appropriate modifications. Modeling can also aid in considering treatments 
and therapies for osteopathies (see e.g. [1]). We feel the type of modeling presented in this work can be 
especially helpful in ruling out ineffective treatments and therapies. 

For future work more sophisticated numerical schemes for solving the level-set equation and computing 
solutions of the model equations will be employed. In particular since remodeling occurs primarily near a 
small region of the interface, a local level-set method such as in [21,28] can be used to speed up computations. 
It may also be appropriate in the future to include further details of the cytokine signaling, something that is 
likely to involve coupling additional differential equations with the existing model. This could call for a need of 
more efficient and accurate numerical algorithms to carry out simulations, particularly important in moving 
to three dimensional problems. While the level-set equation scales well in higher dimensions, numerical 
solutions of other differential equations become much more expensive from a computational standpoint. 
We feel the level-set approach taken here will provide a supplement or alternative to other models of the 
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fundamental process of bone remodeling. Further development in this direction will allow for applications to 
a variety of biomedical issued related to bone remodeling at a variety of scales. 
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